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Abstract. In this paper we propose a numerical method to obtain an ap- 
proximation of Nash equilibria for multi-player non-cooperative games with a 
special structure. We consider the infinite horizon problem in a case which leads 
to a system of Hamilton- Jacobi equations. The numerical method is based on 
the Dynamic Programming Principle for every equation and on a global fixed 
point iteration. We present the numerical solutions of some two-player games 
in one and two dimensions. The paper has an experimental nature, but some 
features and properties of the approximation scheme are discussed. 
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1 Introduction 

The notion of Nash equilibrium is considered one of the most important achieve- 
ments in the last century for its impact on the analysis of economic and social 
phenomena. It is worth to note that the formal definition of Nash equilibria 
[2] has opened new research directions and has attracted the interest of several 
mathematicians to new fields of application. After the pioneering work of Von 
Neumann and Morgenstern (15] , the use of mathematically-based analysis in the 
study of economic sciences has received new impulse from the work of Nash. 

There is extensive literature dealing with Nash equilibria for non-cooperative 
games, however the analysis of Nash equilibria for nonzero-sum multi-player 
differential games is more limited (see e.g. the monographs [3l I13j for an pre- 
sentation of this theory). Moreover, only few papers give a characterization of 
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the value functions for the players in terms of partial differential equations, as it 
was the case for control problems and for zero-sum differential games described 
for example in [5]. More precisely, we know that under suitable regularity as- 
sumptions, if the value functions of a non-cooperative nonzero-sum multi-player 
differential game exist, they satisfy a system of first-order Hamilton- Jacobi equa- 
tions, see [T]. Typically theoretical results about these problems are valid only 
in very special cases and they are essentially limited to games in one dimension 
and with a simple dynamics, see e.g. [H [5J IHl [7] . More important, it is well 
known that the system of Hamilton- Jacobi equations can be ill-posed. To our 
knowledge, there are no theoretical results if the dimension of the problem is 
greater than one. 

From the numerical point of view the situation is even more deceiving since 
few results are available for Nash equilibria in the framework of non-cooperative 
games. In particular, we mention the recent paper [lOj where an approximation 
of (static) Nash equilibria based on Newton methods is presented, and the paper 
[S] where the approximation is obtained via penalty methods. 

Our goal is to construct an approximation scheme for Nash equilibria for non- 
cooperative differential games starting from the characterization obtained via 
Hamilton- Jacobi equations. Following [S], we deal with the system of stationary 
Hamilton- Jacobi equations giving a characterization of the value functions for 
a class of infinite-horizon games with nonlinear costs exponentially discounted 
in time. To this end we will extend to the system the class of dynamic pro- 
gramming approximation schemes studied for zero-sum differential games. The 
interested reader can find in [11] a detailed analysis of the schemes correspond- 
ing to pursuit-evasion games as well as some numerical tests (see also [T2] for the 
origin of these methods and some control applications). At our best knowledge, 
the approximation scheme proposed in this paper is the first for Nash equilibria 
in the framework of differential games. 

The paper is organized as follows. In Section 2 we set up the problem, 
introduce the notations and recall the main results giving the characterization 
of the value functions. In Section 3 we introduce the semi-discrete and fully- 
discrete approximation schemes, and we describe the fixed point iterative scheme 
for the system of Hamilton-Jacobi equations. Some remarks about the fixed 
point algorithm are also discussed. Finally, in Section 4 we present the numerical 
results for some problems in dimension one and two. 



2 Setting the problem 

Let us consider a m-player non-cooperative differential game with controlled 
dynamics 

/ y{t) ^ f{yit),ai{t),...,a„i{t)) , . 

I y{0) = X 

where t > 0, x e W\ / : M" x M''! x ... x M''- M" and Ui : [0, +oo) A, is 
the (open- loop) control associated to the i-th player (z = 1, ...,m) within a set 
of admissible control values Ai C M'''», fcj > 1. We set Ai = {a^ : [0,+oo) — !> 
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Ai , (Xi measurable} for i — l,...,m. In order to simplify the notations, we 
also set a(-) = (q;i(-), . . . , am(0) denote by yxif-', «(•)) the corresponding 

solution of the Cauchy problem ([T]), i.e. the trajectory starting at x which 
evolves according to the strategies a(-) of the m players. 

We consider the infinite horizon problem, where each player has a running cost 
discounted exponentially in time. More precisely, for ? = 1, m, we take Ai > 0, 
■i/'i : M" X M*^' X ... X E'^'" R and we define the cost functional 

J,(a;,a(-)) = / V»(yx(t;a(-)),a(-))e"^'*di, i = l,...,m. (2) 
Jo 

We say that a m-tuple of feedback controls a*{y) = (aj(j/), aj^(y)) (i.e. func- 
tions a* : K" — > Ai, I = l,...,m, depending on the state variable) is a Nash 
non-cooperative equilibrium solution for the game ^ if 

Jj(a;,a*)= min Jj(x, , ...a*_i, a,(-), aJn) : j = l,...,m. (3) 

ai(-)e.Ai 

Note that for every given feedback control : M" ^4^, i = 1, m, and every 
path y(i), i > 0, we can always define the corresponding open loop control 
as ai{y{t)) G Ai, so that, in the above definition, Ji[x,a*) is the i-th cost 
associated to the trajectory yx{t;a*) in which all the optimal Nash strategies 
a* are implemented, namely the solution of 



m^f{yit),a*{y{t))) 
y{Q)^x. 



(4) 



On the other hand the term Ji{x, a^, ...a*_i, ai{-), a*_j_j^, a^) in ([3]) is the i-th 
cost associated to the trajectory yx{t; a*, a*_^, cei{-), a*^■^^, a'^) correspond- 
ing to the solution of 

y{t) = f{y{t),al{y{t)), ...,a*_i,a»(-),a*+i, a™(y(i))) , . 

2/(0) = X 

where only the strategy ai{-) is chosen in Ai- 

The definition of Nash equilibrium means that if the ith player replaces his 
optimal control a* with any other strategy ai(-) G Ai, then his running cost Ji 
increases, assuming that the remaining players keep their own controls frozen. 
The m-tuple a* is then optimal in the sense that no player can do better for 
himself, since he can not cooperate with any other player. 

Let us assume that such a Nash equilibrium a* exists for our game problem. 
For i = 1, ...,TO we define the value function Ui : K" -> R as the minimal cost 
Ji associated to a* . More precisely, for every x € R" and i — 1, m, we set 

Ui{x) = Ji{x, a*) . (6) 

Then, it can be proved that all the Ui satisfy a Dynamic Programming Principle 
and by standard arguments we can derive a system of Hamilton- Jacobi equations 
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for Ml,..., Mm, in which the feedback control a*{x) = (a*(x), a^(a:)) depends 
on the state variable x also through the gradients Vmi(x), \7u„i{x). We get 

XiUi{x) = Hi(x,Vui{x), ...,\7um{x)) a; e R" ,i = 1, ...,m , (7) 

where, for every i — 1, ...,m and for every x, pi, ...TPm G R", the Hamiltonians 
H, : R"+"'" R are given by 

H,{x,pi,...,p^) ^p, ■ f{x,al{x,pi,...,pjn),--,a'^{x,pi,...,p„^))+ (8) 

+ilM{x,al{x,pi, ...,pm), a'^{x,pi, . 
Moreover, for every x G R" and i = 1, ...,m, the following property holds: 

Hi(x,yuiix), ...,VUm{x)) = 

= min \ \7u,{x) ■ f{x,al{x),...,a*_T^{x),a^,a*+:^{x),...,a*^{x)) + (9) 

+ 'ip^{x, al{x), a*_i (x), flj, a*+i(a;), a*^{x)) | . 

We remark that in ^ the minimum is performed among all the control values 
Gi G Ai and not among all the strategies ai{-) G Ai. 



3 Numerical approximation 

This section is devoted to the introduction of the semi-discrete and the fully- 
discrete scheme for the system of Hamilton- Jacobi equations (O , the fixed-point 
algorithm, and to a brief discussion about the algorithm. 



3.1 Semi-discrete and fully-discrete scheme 

Here we propose a numerical scheme to compute the value functions Ui, i = 
l,...,m defined in (jS]). To simplify the presentation and for computational 
purposes, we deal with a two-player game (to = 2), with scalar controls (i.e. 
fci = fc2 = 1) and a dynamics in one or two dimensions {n — 1,2). Moreover, 
we assume the discount rates Ai = A2 = 1. 

In order to discretize the Hamiltonians Hi we use a semi-Lagrangian scheme, 
which usually gives more accurate results than finite difference schemes or oth- 
ers schemes. The reader can find in [12] a comprehensive introduction to this 
subject. As usual we first obtain a semi-discrete scheme, introducing a fictitious 
time step h > 0, which reads as 

ui{x) = min \ -^—ui{x + h f{x,ai,al)) + --^^'ipi{x,ai,a2)\ 

aiGAi + l + h J 

(10) 

M2(x) = min \ ——rU2ix + hf{x,al,a2)) + ——r^2ix,al,a2)> 
026^2 [l + h l + h I 
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where we remind that the control a* — (0^^,02) depends on x and Vui(a;), 
Vu2{x). Note that the time step h can be interpreted as the discrete time cor- 
responding to the approximation of the trajectories starting from x and moving 
according to the dynamics /. 

Now let us consider a subdomain C M" in which we look for the approx- 
imate solutions of the system described above. We discretize by means of a 
uniform grid of size Ax denoted by G = {xi, . . . , cc^v}, where N is total number 
of nodes. For i = 1, 2 we denote by Ui G M.^ the vector containing the values of 
Ui at the grid nodes, i.e. (Ui)j — Ui{xj), j = 1, . . . ,N. Moreover, for i — 1,2, 
let Af be a finite discretization of the set of admissible controls Ai. 

Note that even for x = xj the point z{xj,a*) = xj + h f{xj,a*)) appearing 
in ([TUl) will not coincide with a node of the grid G. Then, an interpolation is 
needed to compute the value of Ui at z (this is the main difference with respect 
to a standard finite difference approximation). First of all, denoting by ||/||oo 
the infinity norm of / with respect to all its variables and choosing 




we guarantee that the point z lies in one of the first neighbouring cells around 
Xj, denoted by I{z). Then, the value Ui{z) is reconstructed by a linear inter- 
polation of the values of Ui at the vertices of the cell I{z) (see [12] for more 
details and [8] for an efficient algorithm in high dimension). Let A(a) denote 
the N X N matrix which defines the interpolation: for j = 1, . . . , A^, the jth row 
contains the weights to be given to the values {Ui)i, . . . , {Ui)N in order to com- 
pute Ui{z{xj,a)). The fully-discrete version of the system (ITU]) can be written 
in compact form as 



where the index j varies from 1 to N . 

We remark that for an actual implementation, the values of U at the boundary 
nodes of the grid must be managed apart, in order to impose boundary condi- 
tions on dVt. In Section |4] dedicated to numerical tests, we will use Dirichlet 
boundary conditions, taking notice of the influence they could have on the nu- 
merical solution inside fi. 

We are ready to describe the algorithm we actually implemented. 

Fixed point algorithm 

1. Choose two tolerances £1 > and £2 > 0. 

Set fc = 0. Choose an initial guess for the values oiUi, i — 1,2 and denote 
them by uf'\ 




(12) 
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2. Forj = l,...,A^ 

(a) Find a* e ylf x A* such that 

' al = argmin (-^ (A{ai,a;) U^"^) 

< 

a*2 = argmin | — f— - ( A (a J, 02) U^^'') 

Note that the search for a* can be done in an exhaustive way, due to 
the fact that the set At^ x is finite. 

(b) If a* is found, go to Step (c), otherwise stop (if more than one Nash 
optimal control is available, we select the first one we find). 

(c) Solve, for i = I, 2, 

(t^^'^'^)^- = YTh (^("*) ^^'0, + (13) 

3. If — t^j^*'''||RN < Si, i = 1,2, then stop, else go to step 2 with 
k <- k + 1. 



1 + h 
h 

1 + h 



3.2 Some remarks about the fixed point algorithm 

We introduce a vector U containing the two value functions, 

U = . . . , {Ui)n, (C/2)i, . . . , {U2)n) e M^A' 

and a vector 5* containing the two cost functions, 

^{a) = {i>i{xi,a), . . . ,^pi{xN,a),ip2ixi,a), . . . ,ip2ixN,a)) e ■ 

We also define a fixed point operator F = {Fi, . . . , F2n) ■ ^"^^ M'^^, given 
component-wise by 

F,{U) = ^{A{a*)U)^ + j^{^{a*)),, j = l,...,27V, (14) 

where, with an abuse of notation, we denoted again by A (a) the block matrix 



A(a) 








A(a) 



In this way the scheme (fT3| can be written as J7 = F{U). 

In the case of a one-player game or two-player zero-sum games with minmax 
equilibrium, it can be easily proven that the corresponding fixed point operator 
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(defined similarly as before) is a contraction map, see for example [TT]. More 
precisely, denoted by G the operator, it can be proved that, for some norm || • ||, 



\\G{U)-G{V)\\<^\\U~V\\, ^U,V. (15) 

In the case of Nash equilibria, the arguments leading to the proof of (jlSp are not 
valid. Nevertheless, we could in principle compute the Jacobian matrix Jp{U) 
of F(J7), i.e.. 

Mm - (2ff>) (16, 

and then study the behaviour of its norm. In this respect, it is important 
to note that the optimal control a* which appears in ([T4| actually depends 
on U , and then A(a*) depends on U . This makes the analytical computation 
of Jf{U) extremely difficult, since we should know the derivative of a* with 
respect to U . Nevertheless, in some particular cases the computation can be 
easy. Indeed, since we assume that A* xA* is finite, it is reasonable to expect 
that small variation of U do not produce a change of a* , i.e. the function a*(U) 
is "piecewise constant". It is worth to note that in the zones of where 
a*{U) is constant, || Jf(C^)||oo is easy computed and we have 

\\JFiU)U^j^^<l (17) 

as in the one-player case or the two-player zero-sum minmax case [llj (the result 
easily follows from the fact that the sum of the elements of each row of A is 1, 
by definition of linear interpolation). Summarizing, it results that F(U) can be 
piecewise contractive in the cxD-norm, according to the following 

Definition 1 A function F : — )> is piecewise contractive in the norm 
II • II if and only if : 

a) F is differentiable a.e. in R^^; 

b) II JF(a;)|| < 1 for all the points x G where F is differentiable. 

From the numerical point of view, we can compute the Jacobian matrix Jp[U) 
of F{U) and check if || Jf(C^)||oo < 1 for some U . Such a computation can give 
some information about the behaviour of F near a tentative solution. Moreover, 
we can plot the graph of some component Fq oi F as a, function of one of its 
variables Ur, again in order to understand the regularity of F. 



4 Numerical tests 

The aim of this section is to provide some numerical simulations in dimensions 
n = 1, 2 in order to give some insights about the well-posedness of system ([7]). 
We believe that such results can be useful to lead future theoretical investiga- 
tions. 
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Let us start by the one dimensional examples analyzed in [S]. The game takes 
place in R with scalar controls ai, 02 S K and it is driven by the dynamics 

/(a;,ai,a2) = ai + 02 . (18) 

The cost functions are 

ipi{x,ai,a2) ^ h{x) + Y ^ « = 1,2, (19) 

where the hi's are smooth functions satisfying |/ii(x)| < C for some constant C. 

For reader's convenience, we recall here a notion of admissible solution for 
the system ^ introduced in [5], for which an existence and uniqueness result 
is available (see [5], Theorem 4). 

Definition 2 A function u — {ui,U2) : R ^ is called an admissible solution 
to the system ^ if the following holds: 

(Al) u is absolutely continuous and satisfies system ^ at a.e. point x gR . 

(A2) u has sublinear growth at infinity, i.e., there exists a positive constant C 
such that, for all x G R, 

\u{x)\ < C(l + . 



(A3) At every point a; e R, the derivative u' admits left and right limits u'{x—), 
w'(a;+). At points where u' is discontinuous, these limits satisfy at least 
one of the conditions 

u[{x-) + u'^ix-) > or u[{x+) + u'2{x+) < . 



Test 1 

We choose the dynamics (|18l) and cost functions with hi{x) = 0, i = 1,2. 
The corresponding system of Hamilton- Jacobi equations is 



(20) 











[ Wl = 


(4- 


f M2 












1 W2 = - 






I u'2 



and admits at least the following three formal solutions: 

Ui{x) = U2{x) = , 

- . ^ / if > 1 - f ^ n 

^^(^) = { -lil-\x\r otherwise, M^) = , 

and ^ 

Ul{x) = ^T,X^ , U2{x) — , 



8 



with optimal controls a* = (— w'l, —u'2), a* = {—u'l, —U2) and a* = (— w'l, —U2) 
respectively. 

It is easy to check that u = (wi,M2) is the only admissible solution among the 
three above. Indeed, the solution u does not satisfy condition (A3) in Definition 
[21 whereas the solution u does not satisfy condition (A2) in Definition [2] 

We solve system (j20p by means of the semi-Lagrangian scheme presented in 
Section IXn Let the numerical domain = [—50, 50] be discretized in 51 nodes. 
Moreover we approximate the sets of admissible controls Ai = A2 = M. with 101 
values uniformly distributed in [—10,10]. For i = 1,2 we set the initial guess 
C/f°'' equal to a large constant. In this case a Nash equilibrium is always found 
and the scheme converges to the admissible solution u = 0. We observe that 
the scheme converges to this solution also starting from the solution u itself or 
from small perturbations of it. Conversely, the two other solutions seem to be 
unstable: even if we choose them as initial guess, the scheme still converges to 
u. We only remark that the case of u requires a larger set of controls, due to 
the large values attained in 51 by the optimal control a* — {—u'l, — M2) — (a;,0). 
Concerning the Dirichlet boundary condition, we have noted that a large con- 
stant does not affect the solution inside the numerical domain and the scheme 
converges quite fast, whereas a negative value leads to the formation of bound- 
ary layers that dramatically slow down (and sometimes hinder) the convergence. 
For example, if we choose the value —10 as Dirichlet boundary condition for both 
Ui and U2, we obtain, after 200 iterations, the results shown in Fig. [TJ 

u, 

2 1 1 1 2 1 1 1 

l[ - 1 - 

CX3XGGGGG0CO00000CCCCCCCCCCCCCCCCCCCCCCCCCG - 0- cffOCCCCCCCCCCCCGCOOOOOOCCOCCCCCCCCCCCCODCO 

o o o o 

-1 - 

o o o o 

-2 h - -2 - 



Figure 1: Test 1. A boundary layer effect for the value functions Ui (left) and 
U2 (right). 

Test 2 

Here we consider again the two-player game described by the dynamics (|18p 
and cost functions (|19|) . where now we choose hi{x) — kix for j = 1,2, with 
constants ki,k2 satisfying fci < < and ki + k2 ^ 0. In this special case the 
solution is explicit, given by 

u{x) = (ui{x), U2{x)^ = (/cix — kik2 — kf/2 , k2X — /ci/c2 ~ ^1/2) • 

Moreover it is the unique admissible solution of system (O, by virtue of The- 
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orem 4 proved in [S]. We choose the same parameters of Test 1 and fci = — 1, 
k2 = 2. Starting with the projection of u on the grid, the algorithm reaches 
convergence immediately. Moreover, we observe convergence also starting from 
small perturbations of u and from a constant initial guess [/(") = C with 

C> max max{ui(x),U2{x)} . 

2;e[-50,50] 

However convergence is not guaranteed for any initial guess. 

Fig. [5] shows the approximate value functions compared with the exact solution 

corresponding to the choice C/'-'^' = C — 150. Following [S], the existence and 




-50 5C -50 50 



Figure 2: Test 2. Comparison between exact solution (solid line) and approxi- 
mate solution (circles), ui and Ui are plotted on the left, U2 and U2 are plotted 
on the right. 

uniqueness result for admissible solutions holds in the more general case in which 
hi{x) is a sufhciently small perturbation of kiX, for z = 1, 2. More precisely, for 
fci, ^2 such that fci < < fc2 and ki + k2 0, it is assumed that 

\h[{x) ^ ki\ < 6 , Ih'^ix) - k2\ < 6 , 

for all a; £ K and a sufficiently small 6 > 0. 

We apply our algorithm also in this case, by choosing S — 2 and 

hi{x) — kix — 5 cos{x) , h2{x) = k2X — 6cos{x) . 

We choose a constant initial guess equal to 100 and the set of admissible controls 
is [—300,300]. The algorithm converges to a reasonable solution shown in Fig. 

13 

Test 3 

Here we present a two-dimensional test with an easy coupled dynamics. We 
choose 

f{x,y,ai,a2) = (02,01) , (21) 
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Figure 3: Test 2. Comparison between exact solution obtained with 5 — (solid 
line) and approximate solution with 6 — 2 (dots and line), ui and Ui are plotted 
on the left, U2 and U2 are plotted on the right. 

with cost functions 

V'^ix,y,ai,a2j - j Q otherwise ' ' ^ ' 

In this game the two players have the same cost function and want to steer the 
dynamics in the unit ball centred in (0,0) where the cost is 0. Considering the 
symmetry of the data, we expect ui = U2- The numerical domain is = [—2, 2]^ 
and it is discretized by 51 x 51 nodes. The sets of admissible controls are 
A\ = A2 — [—1,1], and they are discretized choosing — Af — { — 1,0,1}. 
For i = 1, 2 we set the initial guess U^ ^'' equal to a large constant. Convergence 
is reached in a few hundreds of iterations, the results are shown in Fig. 2] 
The value functions show the expected behaviour, being equal to in the 




Figure 4: Test 3. Approximate solutions Ui (left) and U2 (right). 
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unit ball centred in (0,0) and growing uniformly in every directions. We also 
computed numerically the infinity norm of the Jacobian matrix of F ^ finding 

\\Jf{U)U = ^ 

(see (fT7|) ) in three cases for U: the final solution, the initial guess and an 
intermediate value of the fixed point algorithm. 

Test 4 

The last test is devoted to the investigation of a case where the algorithm does 
not converge. Indeed, choosing 

/(a;,y, 01,02) = (fli +02, ai - a2) , (23) 

■0i(x,y, 01,02) = i = l,2 (24) 

and the other parameters as in Test 3, only some values stabilize, whereas oth- 
ers oscillate. In Fig. [S] we show the surfaces of Ui, U2 obtained after 1000 
iterations. Let us focus our attention on two nodes with different features. 




Figure 5: Test 4. Approximate solutions Ui (left) and U2 (right) after 1000 
iterations. 

In the first node we observe convergence whereas in the second we observe 
an oscillating behavior of the approximating sequence. Let jo G {1,---,2A^} 
be the node corresponding to the central point (0,0) for Ui. We freeze all 
the values oi U = {Ui^U^) except {U)jg, which is replaced by a real vari- 
able s. In Fig. [6] we plot the component Fj^{U) of F{U) as a function of s, 
i.e. J^jo(s) = i^j(,((C/)i,...,(J7)jo_i,s,(J7)jo+i,...,(C/)2Ar). Compared with the 
identity function, it is immediately clear that it is discontinuous and piecewise 
contractive with two fixed points (see Definition [T|) , labelled P and Q (see Fig. 
[5]) . We observe that by our algorithm the value of the node jo converges to the 
fixed point Q. 

Conversely, the value of the node corresponding to the point (0, + Ax) 
does not reach convergence. Similarly as before, we plot the component of the 
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Figure 6: Test 4. First line: the identity function (solid line) and one of the 
component of F as a function of one of its argument (dots and line). The 
function turns to be piecewise contractive with two fixed points P and Q. Second 
line: zoom around P (left) and Q (right). 

function F corresponding to that point, obtaining the result shown in Fig. [71 In 
this case the function is piecewise contractive with no fixed points. The value 
of the node oscillates between two values around the discontinuity, as expected. 
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